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Methods to extract information from the tracking of mobile objects/particles have broad interest in 
biological and physical sciences. Techniques based on simple criteria of proximity in time-consecutive 
snapshots are useful to identify the trajectories of the particles. However, they become problematic 
as the motility and/or the density of the particles increases due to uncertainties on the trajectories 
that particles followed during the images' acquisition time. Here, we report an efHcient method 
for learning parameters of the dynamics of the particles from their positions in time-consecutive 
images. Our algorithm belongs to the class of message-passing algorithms, known in computer 
science, information theory and statistical physics as Belief Propagation (BP). The algorithm is 
distributed, thus allowing parallel implementation suitable for computations on multiple machines 
without significant inter-machine overhead. We test our method on the model example of particle 
tracking in turbulent flows, which is particularly challenging due to the strong transport that those 
flows produce. Our numerical experiments show that the BP algorithm compares in quality with 
exact Markov Chain Monte-Carlo algorithms, yet BP is far superior in speed. We also suggest and 
analyze a random-distance model that provides theoretical justification for BP accuracy. Methods 
developed here systematically formulate the problem of particle tracking and provide fast and reliable 
tools for its extensive range of applications. 

Tracking of mobile objects is widespread in the natural sciences, with numerous apphcations both for living and 
inert "particles" . Trajectories of the particles are to be obtained from successive images, acquired sequentially in time 
at a suitable rate. Examples of living "particles" include birds in flocks [l| and motile cells [1]. Among inert objects, 
nanoparticlcs Q and particles advected by turbulent fluid flow (^-01 provide two important examples. The general 
goal of tracking particles is to extract clues about their dynamics and to make inferences about the laws of motion 
and/or unknown modeling parameters. 

Ideal cases for tracking are those where the density and the mobility of particles is low and the acquisition rate 
of images is high. The non-dimensional parameter governing the stiffness of the problem is the ratio A = Ip^/'^ 
of the typical distance £ traveled by the particles during the time between images and the average inter-particle 
distance \/ p^^'^. Here, p is the number density of particles and d is the space dimensionality. Tracking is rather 
straightforward if A is small: the positions of each particle in two successive images will be relatively far from those of 
all other particles. Trajectories are thus defined without ambiguity. Such a situation is encountered for instances of 
the tracking of nanoparticlcs @. More generally, effective methods are available to identify the assignment (defined 
as a one-to-one mapping of the particles from one image to the next one, i.e., the set of trajectories for all tracked 
particles) that is the most probable [9l-[Tl|. 

The level of difficulty soars as A increases: many sets of trajectories, i.e., many mappings among particles in 
successive images, have comparable likelihoods, see Fig. [1] The dynamics of the particles ought to be described by 
an explicit model, which generally features unknown parameters. The model defines a probability distribution over 
the space of all possible assignments. Contrary to the small A case, the probability distribution is not necessarily 
dominated by a unique assignment. Notwithstanding this uncertainty on the trajectories, it is expected that useful 
information might still be extracted if the number of tracked particles is sufficiently large. The difficulty is that all 
possible assignments must be considered : restricting to the most probable assignment (MPA) generally leads to biased 
inferences (see the sequel). Reliable inference requires summing over all possible assignments with their appropriate 
weights. Problems with large A occur in practice; e.g., even with state-of-the-art cameras, particles in turbulent 
flow and birds in flocks [l| often feature ambiguities in the reconstruction of particles' trajectories. Developing 
systematic methods to tackle these cases constitutes our scope here. 

The plan of the paper is as follows. First, we formulate the problem of particle tracking in terms of a graphical 
model. We then show that an exact and rapid algorithm for summing over all possible assignments is unlikely to 
become available, as such an algorithm could equivalently compute the permanent of a non-negative matrix, a problem 
that is well-known to be #P-complete [l^. An approximate message-passing Belief Propagation (BP) algorithm 
[l6j is then introduced, employed and tested. We also introduce a simplified model where analytical results and 
a quantitative sense of the BP approximations are obtained. The Results section presents numerical simulations 
comparing results of BP and the inference based on the MPA. For the sake of concreteness, we consider the case of 
particles passively transported by a turbulent flow, but the methods arc quite general and can be applied to other 
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FIG. 1: A concrete example of particle tracking, with N = 400 particles moving from their original positions (red circles) to 
new ones (blue diamonds). In the left figure, two consecutive images are superimposed to facilitate comparison of the successive 
positions. Particles are transported by a turbulent fluid flow with local stretching, shear, vorticity and diffusivity parameters 
a* = 0.28, 6* — 0.54, c* — 0.24 and k* = f .05 (see eq. (O). We focus on turbulent transport because of the challenges it poses, 
yet the methods we develop are quite general. The right figure shows the actual motion of each particle. Evidently, the simple 
criterion of particle proximity fails to pick the actual trajectories, and the mapping of the particles between the two images 
is intrinsically uncertain. Nevertheless, the inference algorithms described here rapidly yield excellent predictions (a = 0.32, 
b = 0.55, c — 0.19 and n — 1.00) for the parameters of the fiow. 



situations as well. 



I. MODELS 



Tracking of particles as a graphical model. Graphical models provide a framework for inference and learning 
problems widespread in machine learning, bioinformatics, statistical physics, combinatorial optimization and error- 
correction [3, uB, The assignment problem involved in the tracking of particles is conveniently recast as a 
weighted complete bipartite graph (see Fig. 2). Nodes are associated with the N particles in each of two successive 
images, their positions being denoted Xi and , respectively. These 2N position vectors («,j = 1, . . . ,N) constitute 
the experimental data provided. We suppose that a model for the dynamics of the particles is available and features 
a set of unknown parameters 6. Edges between nodes of the bipartite graph are weighted according to the likelihood 
that, according to this model, a particle moves from the initial position Xi to the final position . Specifically, the 
formula for the likelihood of an assignment among (non-interacting) particles in two images reads as follows : 

£(W,0) = c({a}) n \pi{x.,y^e)\^^ . (1) 

(ij) 

The Boolean variable af indicates whether the particles i and j are matched {af = 1) or not (cr^ = 0). The set of 
the N'^ variables aj is denoted by {a}. The constraint function C ({cr}) = ^ ^{^ l) Y\i ^ ' involving 

Kronecker 5 functions, enforces the conditions for a perfect matching, i.e., a one-to-one correspondence between the 
particles in the two images. (Situations where the number of particles in the two images can differ and/or positions of 
the particles are uncertain are accommodated within the same formalism presented in the sequel; see Supplementary 
fnformation (SI).) The quantity P/ is the transition probability that a particle at position Xi travels to yj in the time 
A between images. The transition probabilities carry all of the information about the model for the dynamics of the 
particles. 

The tracking inference problem that we address here is to provide fast and reliable estimates of the model's unknown 
parameters 6, which enter the likelihood of the trajectories via eq. ([J). 

The simplest method to infer the unknown parameters is first to identify for any the MPA, i.e., a configuration 
{cr} satisfying the constraint C({cr}) — 1 and having the highest likelihood ([T]), and then to maximize the resulting 
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FIG. 2: The complete bipartite graph for the tracking problem. Nodes (squares) denote particles in two consecutive images. 
Edges carry weights P^ (the likelihood that a particle travels from the initial position Xi to the final position y-' ) and Boolean 
variables af (indicating whether the nodes i and j in the two images correspond to each other (cr = 1) or not {a = 0)). Conflicts 
arise from the constraint that a valid assignment is a one-to-one mapping of the particles in the two images (expressed by the 
constraint C{{a}) in eq. ([T])). 



likelihood with res£ect to 9. Exact polynomial algorithms [T^] solve the first MPA task. The surprising remark 
recently made in [ll|, [l^ is that the auction exact algorithm can be reformulated as a message-passing Belief Prop- 
agation (BP) scheme. (Surprise stems from the fact that BP usually differ from the exact solution of a problem if 
loops are present in its graph.) 

It is however expected — and confirmed by results described shortly — that the MPA provides reliable inferences only 
for small enough values of the stiffness parameter A defined in the Introduction. As A increases, assignment-dependent 
entropic factors become important, and MPA inferences deviate from the actual values of the parameters 9. This 
deviation is understandable because the Bayesian probability distribution for the parameters, assuming uniform prior 
probability for the assignments, involves the full likelihood ([T]) marginalized over all possible assignments : 

P(0|a;„y^)cx^£({a}|0)EEZ(0) . (2) 

As the stiffness parameter A increases, many assignments become likely and the sum ([2]) is not dominated by the 
MPA. 

Summing over all possible trajectories. In the vast majority of cases, no exact algorithm is available to sum 
over all possible states of a graphical model. Equation is no exception, as it can also be seen as a sum over 
all possible permutations of the lower indices i into the upper indices j. It is then recognized that computing the 
likelihood Z (also known as partition function in statistical physics) is equivalent to computing the permanent of 
the matrix , a well-known #P-complete problem [Tsj . However, the matrix in our sum Q is non- negative 
and the permanent of non-negative matrices was the first #P-complete problem discovered to be solvable by a Fully 
Polynomial Randomized Approximation Scheme (FPRAS) The complexity of the original FPRAS algorithm is 
0(iV"). We significantly accelerated (to 0{N^)) and simplified the original Markov Chain Monte-Carlo (MCMC) 



algorithm of [17| without observable deterioration of its quality (see SI). In the Results section, we use this simplified 
version to assess the accuracy of our BP approximation while in the SI we compare performance of BP to the MCMC 
scheme. 

Among the possible approximations to compute the permanent of a matrix, Belief Propagation (BP) has a special 
status because of its aforementioned exactness for the MPA problem [ll|. Moreover the BP algorithm is very fast, 
scaling as 0{N^) in its basic form and linearly if, for each particle, only a limited number of nearby particles is 
considered. We shall therefore pursue the development of BP to approximate the sum in ^ and then assess its 
validity through numerical simulations and the simplified "random distance" model discussed shortly. 

The starting point of the BP approach is the remark that the convex KuUback-Leibler functional 

^{5(W)}^^6(W)lnM^, (3) 

{a} '-KX^i) 
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has a unique minimum at b{{a}) = C{{a})/Z (under the normahzation condition X]{cr} ^({"'D " 1)' where C is 
defined in eq. ([1]) and Z in eq. ([2]). The eorresponding value of the functional T is the log-likelihood (free energy), 
J- = — In Z. This remark constitutes the basis for variational methods (see [l^), where the minimum of the functional 
is sought in a restricted class of functions. BP (and the corresponding approximation for the free energy, named in 
[l6j after Hans Bethe) involves an ansatz of the form 

where each vector <t,; = {af |j = 1, • • • , N} can be any of the N possible vectors (0, • • • , 0, 1, 0, • • • ,0) having exactly 
one nonzero entry, and cr-' is defined analogously. This ansatz is motivated by the fact that the probability distribution 
on a graph with a tree structure, i.e., without loops, takes exactly the form (|4]). The quantities bi, V , and 6^ are 
called "beliefs". In the absence of loops, represents the probability distribution for the single Boolean variable cr^, 
and bi represents the joint probability distribution of the components of cr^, which appear in a constraint. Without 
loops, the beliefs automatically satisfy the conditions of marginalization, viz., 6^(cr^) = X^o- \[t^ bi{cri)^ where the sum 

is performed over all possible values of the vector cr,; having the specified value of its component cr^ . 

In the presence of loops, the expression ^ is no longer exact, and the marginalization conditions must be imposed 
as additional constraints. The important point demonstrated in [l^ is that minimizing the functional ^ for the class 
of functions (j4|) (under the marginalization and the normalization conditions) yields the message-passing formulation 
of BP, which was heuristically introduced by Gallager for decoding of sparse codes [l^ . 

Most commonly, the BP equations are formulated as an iterative message-passing scheme [H,[23|, where the message 

jj^^o (respectively, h^'') is sent from from particle i in the first (second) image to particle j in the second (first) 
image. In the case of the matching problem, the messages are determined by solving the following equations: 

t^' = -iln^/f e^^^-' ; h^^^ = -^In^^e^^^^^ . (5) 

The "inverse temperature" /3 can be set to unity, but it is usefully retained to show that the limit /? — >■ oo yields the 
exact solution for the MPA problem [Tl| . 

To solve the BP equations ([5]), the messages are randomly initialized and iteratively updated in order 

to find a fixed point of the message-passing equations ([S]). The Bethe free energy then reads 

The Bethe free energy J- bp, evaluated at the fixed point of the BP equations ([S]), provides an estimate of the exact 
log-likelihood — In Z{6) defined through ([2]). Because the most likely set of parameters 9 minimizes — In Z{9), we 
seek parameters that minimize the estimated free energy. We perform this minimization using Newton's method in 
combination with message-passing: after each Newton step, we update the messages [{h}, {/i}) according to ^ and 
the current set of parameters 6. As this combined update takes iV^ steps and the number of combined step which 
led to convergence is A^-independent, the running time scales as 0(N'^). Note that the running time can be further 
reduced to 0{N) neglecting the contribution of edges with very small probability P/. i.e. diluting the fully connected 
bipartite graph. 

A simplified model to understand the BP approximation. The BP approximation is exact only if the un- 
derlying graph is a tree, which is not the case for our fully connected bipartite graph in Fig. 2. In the Results, we 
empirically assess the validity of the BP approximation through numerical simulations. To supplement the numerical 
evidence, we introduce here a simplified "random distance" model for which analytical results can be obtained and 
used to understand the nature of the BP approximation. 

The "random distance" model is defined as follows. First, we decouple the A^^ distances between particles i 
and j by assuming that they are independent among each other. We then assume that one permutation tt* of the 
lower indices i into the upper indices j has a special status, while all other distances are drawn independently 

at random from a given distribution. Namely, the A^ distances dl~^^ 's are distributed as a Gaussian (restricted to 
positive values) with variance k* = 0(1). For each of the other N{N — 1) pairs («, j), the distances dj are independent 
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FIG. 3: A realization of a 2-dimensional flow with a* — b* — c* = 1, k* — 0.5 and A'' = 200 particles. Left: The BP Bethe 
free energy as a function of the diffusivity k and the shear b, where every point is obtained by minimizing with respect to the 
stretching a and the vorticity c of the flow. Right: The same free energy in a contour plot, showing the maximum close to 
b — 1 and k — 0.5. The maximum is achieved for asp = 1.148(1) 6sp — 1.026(1) cbp = 0.945(1), kbp ~ 0.509(1), where the 
parentheses indicates numerical error on the third digit. 

random variables drawn uniformly in the interval (0, A^). Units of length are chosen so that the typical interparticle 
distance is set to unity. Note that any distribution of dj for j ^ n* with a vanishing derivative at the origin would 
give the same solution. Indeed, the crucial property is that, for each particle i, the number of distances that are 
comparable with the diffusion length scale y/Ti* is 0(1). 

The interest of a special permutation tt* is that we can inquire about : learning k* if tt* is supposed unknown; the 
relevance of entropic factors for the partition function; and the status of the BP approximation. The same questions 
arise for ^ in the original problem. Note that if there were no special permutation, then we would obtain the 
random-link model considered in |22h24| . Our "random distance" model can be solved exactly in the thermodynamic 
limit using the replica method, as in |22l. [23j . or using the cavity method, as in [2^. The main result (see SI) is that 
the BP expression of the free energy for the "random distance" model is exact in the limit N ^ oo, despite of the 
short loops in the graph of the model. The argument to prove this result goes as follows : (a) the contribution to the 
partition function ^ from those permutations of the j indices that contain distances larger than 0(1) is negligible; 

(b) as each node has only a fraction 0(1) of its A^ distances being 0(1), the underlying graph is effectively sparse; 

(c) because sparse graphs arc locally tree-like and correlations in the matching problem decay very fast on trees, it 
follows that the BP approximation is exact in the thermodynamic limit. We have formalized these statements within 
the replica and cavity methods (and rigorous local weak convergence methods are probably also applicable, as for the 
random link model |25|). 

The asymptotic exactness found in the "random distance" model means that errors made by BP are caused by 
correlations among the inter-particle distances. In a smooth flow (see the sequel) particles close to each other in the 
first image will also be near each other in the second image; moreover, the four distances among the two positions 
in the first image and the two positions in the second will also be small, or more generally correlated. This effect is 
more important in lower dimensions. Indeed, BP inferences turn out to be better in the three-dimensional (3D) case 
than in 2D and rather inaccurate in ID (with maximum relative error about 60%). 

Our replica calculations also show that the "random distance" model presents an interesting phase transition at the 
diffusivity Kc ~ 0.174. For k* < Kc, the MPA ttmpa is identical to the special one tt* with high probability, whereas 
for K* > Kc the overlap (defined via the Hamming distance) between the most likely assignment ttm pa and the special 
one TT* is extensive, i.e., 0{N). The comparison with the finite-dimensional case is discussed in the Results section. 
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FIG. 4: The BP estimate of the log-hkelihood —J- bp and the Monte-Carlo Markov-Chain estimate plotted vs. the diffusivity 
K with K* — 0.1, 1, 10 for = 100 particles diffusing in 3D. Although the BP log-likelihood is significantly lower than the 
MCMC one, the estimate of the diffusivity is extremely good. The vertical lines mark the diffusivity computed based on the 
knowledge of the actual displacements. The BP algorithm compares favorably in estimation of the maximum with the basically 
exact Markov Chain Monte-Carlo algorithm, while being far superior in speed (for details on speed comparison see SI). 



II. RESULTS 



Our analysis has been quite general so far. To concretely assess the validity of BP, it is now necessary to specify the 
model appearing in the likelihood ([1]), namely the probability P/ for the transition from position Xi to in the image 
acquisition time A. We decided to focus on the tracking of particles in turbulent flow for three reasons. First, the 
problem is highly relevant as an important part of modern experiments in fluid dynamics is based on the tracking of 
particles either in simple [^-Q] oi complex (26| flows. Second, all the algorithms used so far for reconstructing the flow 
from images of multiple particles have ignored the probabilistic structure of the possible assignments. As discussed in 
the review [l^ . the general approach is to search for a single matching. Criteria based on proximity and/or minimal 
acceleration identify for each particle its "best" mapping in the successive time-image. Conflicts where two or more 
particles in the first image are assigned to the same particle in the second image are resolved by various heuristics. 
The simplest option is to give up on those situations where a conflict arises; the most elaborate solution is to compute 
the assignment with the minimal cost in terms of proximity or acceleration. The bottomline is that one is always 
left with a single assignment, which leads to predictions that are effective at low density of the particles but rapidly 
degrade as their density increases [l^ . Third, the laws of motion of the particles arc well-known. Indeed, if particles 
are sufficiently small and chosen of appropriate (mass) density, their effect on the flow is negligible and they are 
transported almost passively. The (number) density of particles is usually rather high and a single snapshot contains 
a large number thereof. The reason is that the smallest scales of the flow ought to be resolved. Furthermore, turbulence 
is quite effective in rapidly transporting particles so that the acquisition time between consecutive snapshots should 
be kept small. Modern cameras have impressive resolutions, on the order of tens of thousands frames per second. 
The flow of information is huge: ^ Gigabit/s to monitor a two-dimensional slice of a (IOcto)^ experimental cell with 
a pixel size of 0.1mm and exposure time of 1ms. This high rate makes it impossible to process data on the fly unless 
efficient algorithms are developed. 

Tracking of particles in a turbulent flovif. The likelihood P/ in eq. ([T|) for the transport of particles in turbulent 
flow is obtained as follows. Tracked particles are supposed to be at reciprocal distances smaller than the viscous scale of 
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FIG. 5: Scatter plot of the diffusivity estimated by BP and by MPA vs. the actual value of the diffusivity. Diffusion takes 
place in 3D, with displacements generated using n* — 1. The actual value of the diffusivity factual is computed from the actual 
trajectories and is subject to statistical fluctuations. The number of tracked particles is A'^ = 100 (red) and N — 400 (blue) 
using 1000 measurements. The BP predictions correspond to the maximum of the log-likelihood, as approximated by the Bethe 
free energy ((6]) discussed in the text. Note the strong underestimation of the MPA estimate, to be contrasted with the cloud 
of BP predictions centered around the correct value of k* . 



the flow. It follows (see [21| for a review) that the position ri(t) of the i-th particle evolves according to the Lagrangian 
stochastic equations : = s-ri+^i. Positions are measured with respect to a reference point and s is the tensor of the 
velocity derivatives. In two-dimensional incompressible flow a = s^x — '"Syy is the rate of stretching, b ~ {sxy + Syx) /I 
is the shear, and c = {sxy — Syx)/2 is the vorticity. The stochastic term ^i(t) is the zero-mean Gaussian Langevin 
noise, describing molecular diffusivity, defined by its correlation function: {{^i)a(ti){^j)p{t2)) = 2KSijSapS(ti — t2). 
The Greek indices refer to space components. The transition probability corresponding to the previous transport 
process is Gaussian: 

P^{x„y^) = (deti\/)-5exp(-ir"(Af-i)"'3r'3); (7) 
r = yJ" - W{A) X, ; (8) 
M = KW{A)\S^W-^{t)W-^''^{t)At] W^^(A), (9) 

where W(t) = exp(i s) (35| and A denotes the image acquisition time. The general problem of estimating the unknown 
parameters 9 now takes the special form of inferring the components of the tensor s and the diffusivity k. 

We tested the BP inference method on numerically simulated data. To compare different system sizes, we placed 
the N particles at random in a d-dinicnsional box of size L = N^/'^, i.e., the average density equals unity. Particles are 
then displaced independently following the probabilistic distribution ([7]) with a set of parameters a*-/V~^/'^, b*N~^^'^, 
c*iV-i/'i^ K*. The timescale was chosen as the acquisition time A = 1 and the parameters a, b, c of the flow were 
rescaled by N~^^'^, so that the particle displacements in the acquisition time are 0(1) for all choices of TV. Fig. 3 shows 
the BP free energy ([H]) as a function of the shear b and the diffusivity k ioi N = 200 particles. The curvature around 
the minimum of the exact free energy is inversely related to the statistical error in the estimation of the parameters. 
Figure 3 clearly shows that the most problematic parameter is the diffusivity k, as confirmed by all of the numerical 
simulations we performed. 

In the next three figures, we concentrate on the purely diffusive regime (where a* — b* = c* = 0); subsequently, we 
return to the general case. Note that light scattering experiments provide an established measurement method for 
biological systems j27l - l30j . The technique is not commonly used in fluid dynamics experiments because the dispersion 
of the particles is much faster and the typical illumination level is too weak. 

In Fig. 4 we compare the BP results with results from the fully polynomial randomized approximation scheme 
based on the Monte Carlo Markov Chain (MCMC) method. As computed by MCMC (which is guaranteed to be a 
FPRAS), the maximum of Z(k) coincides with the true value k*. From the large deviation theory interpretation of 
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it follows that the statistical error of the estimation of k is cr = \I\^NT"{k)\. Figure 4 confirms the expectation 
that the curvature at the minimum of the free energy decreases with the diffusivity constant k. 

In Fig. 5 we compare the estimates of the diffusivity using BP (kbp) and using MPA (km pa)- The actual value 
^actual (respectively, the MPA estimate km pa) is computed as the mean-square displacement of the particles on the 
actual (resp., the most probable) trajectories of the particles. The key conclusion to be drawn from Fig. 5 is that 
MPA largely underestimates the diffusivity, whereas the BP method is accurate. 

Fig. 6 gives a quantitative sense of the BP accuracy as a 

function of the diffusivity. The upshot of the curves is that MPA gives accurate estimates only at extremely low 
diffusivitics. Conversely, as the diffusivity increases and the overlap among possible assignments becomes important, 
the quality of MPA predictions degrades very rapidly. Results from a similar study for 2D flows is presented in Fig. 7; 
vectorial parameters are again found to be computed efficiently by the BP method. 
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FIG. 6: The relative error ^/'^{k.bp/mpa ~ KactuaO^Z-ft^/^actuai in the estimates of the diffusivity over K measurements vs. 
the actual value of the diffusivity k*. Circles and squares refer to Belief Propagation (BP) while triangles and lozenges refer 
to the Most Probable Assignment (MPA). fCactuai is the actual value of the mean-square displacement of the particles, i.e., it 
includes fluctuations around k* due to the finite number A'^ of particles. The data are averaged over 1000 (for = 100) and 
250 (for N — 400) realizations and compared to relative the statistical error \j2ldN (dashed horizontal lines). The case of 
two-dimensional (2D) diffusion is shown on the left side, and the 3D case is on the right side. The top figures are zooms into 
the low diffusivity region. 



Finally, Fig. 6 (top panel) indicates that the phase transition in the exactness of the optimal assignment, which was 
previously found for our "random distance" model, is smeared out in the finite dimensional case (or it happens only 
at K < 0.01). This effect is traced to the fact that, in finite dimension, there are always several pairs of particles at 
distance o(l) that get confused with the diffusion (when k* = 0(1)). It follows that the probability for ttmpa — tt* 
is always smaller than unity. 
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FIG. 7: Scatter plot of the parameter estimations using the BP method in the case of a 2D incompressible flow with the rate 
of stretching o* = 0.1, the shear b* = 0.5, the vorticity c* — 0.2 and the diffusivity k* — 1. Red (blue) points refer to the case 
of a number of tracked particles A'^ — 100 (A'' = 400) and the number of measurements is 50. 

III. DISCUSSION 

The message-passing algorithms discussed here were shown to ensure efficient, distributed and accurate learning of 
the parameters governing the stochastic map between two consecutive images recording the positions of many identical 
particles. The general method was illustrated in a model relevant for the tracking of particles in fluid dynamics. It 
was shown that parameters of the flow transporting the particles could be efficiently and reliably predicted even in 
situations where a strong uncertainty in the particles' trajectories is present. 

We introduced and compared two techniques to approximate the likelihood that the dynamics of the particles is 
compatible with the displacements observed in the experimental snapshots. The first is based on finding the most 
probable trajectories of the particles between the times of the two images. The second corresponds to evaluating the 
probabilistically weighted sum over all possible trajectories. The latter is a #P-complete problem and its solution 
is approximated by Belief Propagation (BP), as implemented via a message-passing algorithm. BP was shown to 
become exact for the simplified "random-distance" model we introduced here. In general, the effect of loops in 
the graphical model for the tracking problem remains nonzero even in the thermodynamic limit of a large number of 
tracked particles. Preliminary analysis of the loop corrections to BP did not display any immediately visible structure, 
yet detailed analysis of this point is left for future work. Another interesting direction is the development of learning 
algorithms (both MCMC and message-passing) specifically designed to provide estimations of appropriate observables, 
e.g. the sum of the square of the distances traveled by the particles, whence the parameters of the dynamics, e.g. the 
diffusivity, can be estimated. This could lead to further reductions in the computational time and it will be of interest 
to test whether the superiority in speed we found here for BP as compared to the FPRAS Monte-Carlo scheme (see SI) 
still holds. The price of a single observable is that the log-likelihood curves in Fig. 4 offer more complete information, 
namely a systematic way to gauge the error bars on the inferred parameters. 

The algorithms presented here can be carried over to the tracking of other types of "particles", e.g., those of 
biological interest. Motile bacteria in colonies, eukaryotic cells or fluorescent bio-molecules provide relevant examples. 
As it was stressed here, our methods are poised to deal with dense conditions where a strong overlap among the 
various particles' trajectories are present. An additional use of the techniques introduced here is to compare different 
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models of transport, e.g., purely diffusive, directed, active, etc. The Bethc expression given by the BP approximation 
can be taken as an approximation for the log-likelihood of the various models, which are then compared by standard 
model selection tests. The validity of the various models postulated for the dynamics of the tracked objects can thus 
be quantitatively compared. 

Our main conclusion is that the BP method gives accurate results and its computational burden is comparable to 
identifying the most probable trajectories. The accuracy of BP was shown to compare extremely well with exact results 
and to improve rapidly as the dimensionality of the problem increases. The BP-based technique allows generalization 
to reconstruction of a multi-scale from particle images in two sequential snapshots. BP can also be adapted to 
the case where trajectories of the particles are reconstructed from their positions in several (> 2) images (see SI). 
In conclusion, the formulation of particle tracking as an inference problem permits tackling it systematically and 
introducing message-passing methods that are highly effective in diverse applications. 
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Appendix A: Implementation Details 

Generating the data To generate the set of positions Xj in the original image we proceed as follows. We place 
N particles with positions selected uniformly at random in a d-dimcnsional box of size N-^^'^, i.e. each of the d 
coordinates of any of the N particles is i.i.d. on the interval {O^N-^^'^). We choose the set of governing parameters, 
6 = {k,s), describing the diffusivity and the flow transporting the particles, in such a way that the displacements 
of the particles in the acquisition time of the images (chosen as the unit of time) is of order unity (we remind that 
the unit of length is chosen as the typical inter-particle distance) . Namely, since the typical displacement will be of 
the order sr and r = 0(iV^/'^), we need that each component of the matrix s be 0(A^^^/''). We therefore performed 
the rescaling s = s*N~^^'^, with s* = 0(1). To produce a stochastic map for the N particles from the original 
image to their respective positions in the successive image, we generate a d-uple of coupled Gaussian variables with 
covariance matrix according to eq. (8-9) and then output yi = ri + Wxi, where the evolution operator W is 

defined in eq. (9). Note that this setting guarantees that, for sufficiently large N , the number of particles leaving the 
box is much smaller than those staying in the box. For the purpose of our validation, particles which leave the box 
are treated as if the box did not exist. 

Computing the actual estimate, ©actual Given the set of positions Xi and and the actual permutation tt* 
(in practice, tt* can be taken as the identity without any loss of generality), we consider the likelihood eq. (1) to be 
a function of the set of parameters 6. The point in the 6 space where the likelihood achieves its maximum is called 
"the actual estimate" and it is denoted as ©actual- The estimate ©actual is found using the Newton's method. Note 
that finite size effects make that ©actual typically differ from the respective original value, ©*, by deviations that are 
0(Ar-i/2). 

Computing the MPA estimate, ©mpa Computation of the Most Probable Assignment (MPA) estimate is split 
into two parts. First, for the given set of parameters © we look for the permutation, ttmpa, which maximizes the 
likelihood eq. (1); second, we maximize the resulting likelihood with respect to © using the Newton's method. Exact 
algorithms are available to achieve the first task j9l4ll| , and we utilize in our simulations the zero-temperature Belief 
Propagation algorithm of [Tl| . as the latter is considerably faster than other alternatives. In the zero temperature 
limit, /3 — >■ oo, the BP equations (5) become 

V^' = - max(/i'=^* + In P*^) , = - m&yi{h^~'^ + In P/ ) . (Al) 

The MPA matching is reconstructed from the fixed point of cq (jAip as follows 

7ri(i) = -argmaxj.(/i'''^' + InPf ) , 7r2(fc) = -argmaXj(7i'^'' + Ini^'') . (A2) 

We initialize the messages {/i}, {/i} at random, and solve eqs. (|A1[) iteratively till the consistency condition 7r2 [tti (i)] = i 
for alH = 1, . . . , iV is satisfied. Note that this scheme is faster than enforcing the convergence of the BP equations 
(jAip to a fixed point. A •priori we could reach a configuration of messages where the condition is satisfied and the 
corresponding configuration is not the best matching because the fixed point was not reached yet. However, in reality 
we never observed this situation when the number of particles is sufficiently large. 

Computing the BP estimate, ©bp The task here is to find the minimum of the Bethe free energy eq. (6) 
over the set ©, where solve eq. (5) at /3 = 1 and given ©. Our implementation of this task is as follows. 

We initialize the algorithm with randomly generated {/i}, {/i}, solve BP equations iteratively and then minimize the 
resulting Bethe free energy using the Newton's method. Each time the Newton's method calls for the value of the 
free energy at a certain value of parameters ©, each BP message is updated on average m times and then the value of 
the Bethe free energy is computed. We typically used m = 5 or to = 10. In general this small number of iterations is 
not sufficient to reach the fixed point for each value of the parameters, however, we found empirically that reaching 
the fixed point at every step is unnecessary as long as the Newton's method finally converges (which was always the 
case in our implementation). 

Mixing the BP iterations with the Newton's method updates accelerates the algorithm significantly. The running 
time of this mixed implementation is quadratic in the number of particles. Moreover, the following procedure allows 
to attain a linear scaling without much of a quality sacrifice: for each particle i, we retain only those j's such that P/ 
exceeds a predefined small constant C . As the size of the system grows, this truncation identifies only 0(1) number 
of possible partners for each particle resulting in the linear scaling for the algorithm running time. 



Appendix B: Cavity solution for the random distance model 

We give here a brief exposition to the cavity solution of the "random distance" model defined in the main text. 
More details about the method itself, applied to a different but related model of random matching, can be found in 
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In the random distance model the hkchhood (1) becomes 



Ci{a}\K) = C{{a}) n \Pm\^ 



1 


■ 1 " 







(Bl) 



where P/(d^|K) = ^i=e~^^(''i)^ Note that X^r, ~ = ^ 1) and thus the term with 1/A^ in eq. ([BT]) is 

contributing as a multipUcative constant and can be ignored. Hence, the behef propagation equations (5) hold for the 
random distance model, the Bethe free energy is computed according to eq. (6), and the respective expression for the 
ground state energy is 

E^y^ max(0, h'^^ + h?^' + In P/ ) - V max(/i^^' + In P/ ) - V max(7^'^^ + In Pf ) . (B2) 

The belief propagation solution applied to an instance of the random distance problem becomes asymptotically and 
statistically exact in the thermodynamic limit — > oo. This statement is equivalent to exactness of the so-called 
replica symmetric solution made in (23 - 123 . [sH ] for a different but related random matching model. We tried to detect 
instabilities towards replica symmetry breaking in the random distance model, but failed to find any. This empirical 
evidence leads us to conjecture that one may be able to extend the rigorous local weak convergence method of [25| or 
independent combinatorial proof of [32| to the random distance model. 

Population dynamics solution Population dynamics [s^ is a general method to solve the belief propagation 
equations on a infinite sample without actually generating or storing the whole sample. The assumption here is that 
all the quantities of interest are self-averaging, i.e., their values are asymptotically equal to respective averages over 
large sample. For the matching problem, a trick introduced in [2^ . [sij exploits the effective sparseness of the fully 
connected factor graph. The main observation is that the best matching has energy proportional to N whereas a 
random permutation has energy proportional to N'^ , and only permutations with energy proportional to TV contribute 
to the log-partition function and other quantities of interest. Edges with 0{N) weights almost never contribute to 
these configurations. The number k of edges adjacent to a variable whose weights are smaller than a certain constant 
c is described by the following Poissonian distribution with mean c: 

Additionally to these k edges we always consider an edge connecting the variable under consideration to its actual 
image. We thus end up with a bipartite factor graph with the connectivity degree 1 + /c, where k fluctuates according 
to eq. (|B3p . As c — > oo we recover the original problem, but in practice moderate values of c are sufficient to achieve 
numerically the asymptotic regime. 

Let us now turn to the cavity equations explaining the relations among the probabilities of the messages h^^-' , h 
at the fixed point of eq. (5) on a very large graph. In fact, there are two kinds of BP messages, these describing 
transitions between i and its actual image and also transitions between i and nodes which do not correspond to i's 
actual image. Our notation for the two distributions are Q{h) and Q{h), respectively. The resulting cavity equations 
are 

Y[ dd, U{d,) / n <^h^ Q{K) 5{h - T{{hi\)) , (B4) 

K i—1 i—1 

/A, rt 
ddo p(do) n d'^' ^('^') / d'^o Q^^^) n Q^'^^) - •^({'*^}' '^o)) ■ (B5) 

K i=l i=l 

The quantity J^{{hi}) satisfies eq. (5). Let us recall that P{d) ~ —^==e^2h''^ is the distribution of displacements; 
and U{d) — 1/c ioi < d < c, and U{c) = otherwise. 

Population dynamics is the method we utilized to solve the cavity equations (jB4IIB5[) . In this method both dis- 
tributions Q{h) and Q{h) are represented by a pool of A^pop numbers initialized at random. In one sweep on the 
population dynamics we repeat iVpop times the following procedure: 



• Iteration of eq. (jB4|): Draw random number k according to eq. (|B3p . draw k random numbers d,; uniformly 
at random from interval (0, c), and draw k random numbers from the pool representing the distribution Q{h). 
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Then use eq. (5) to compute a new h and substitute it on a place of a random element in the pool representing 
Q{h). Note that in the case when k = from the definition of message h, the value of h has to be large h oo. 
In practice, we take h = h„ 



'-max ■ 



• Iteration of eq. (jBSj: Draw random number k according to cq. (|B3p . draw k random numbers di uniformly 
at random from interval (0, c), and draw k random numbers from the pool representing the distribution Q{h). 
Draw one random number d from the Gaussian distribution with variance k* , and one random number from the 
pool representing Q{h). Then use eq. (5) to compute a new h and use it to replace a random element in the 
pool representing Q{h). 

We iterate the procedure till convergence. Knowing the convergent Q{h) and Q{h) allows us to evaluate various 
sample averaged quantities of interest, such as the density of the free energy (6). Note that to calculate the first term 
in cq. (6), the term needs to be split in two parts, correspondent to the actual {ij) pairs (particle and its image) and 
"confused" pairs. In practice we compute each of the terms in eq. (6) iVpop times, each correspondent to new choice 
of randomness. Then, we do few (around 10 in practice) sweeps to update the pools representing Q{h) and Q{h) and 
repeat the measurement to eliminate sampling errors. Similar scheme also applies to evaluating average of the second 
derivative of the free energy. 

Let us also mention for completeness that the system of cavity equations (jB4IIB5[) allows an elegant and compact 
direct representation in the asymptotic c — > oo limit, in the spirit of (22l. [23j where this trick was originally proposed 
for the random matching model. Following [23| we define a generalized Laplace transform as 

/•oo 



AhQ{h)e-^ , e'^^'^ = / <lhQ{h)e-^ . (B6) 

} J —oo 

Applying this transform to eq. (jB4IIB5P we arrive at the following set of closed equations for G{1), G{1): 

^y^'""^"^ E „ .„ e^^'+"^gp(/3) : (B7) 

^p!(p-l)! 



G{1) = 0(0 -log 

where 



/oo _ jfo , -|\T|-1 



(B8) 



Note, that cq. (|B7[) with G ~ G was originally derived in [22|, [2^ using the replica method for the random matching 
problem, while the log-correction on the right hand side of eq. (jB8[) is specific to our problem, representing bias 
brought into the problem by the existence of the special permutation tt* . 

Results Let us briefly discuss the results obtained by the population dynamics numerical evaluation of the cavity 
equations cq. (|B4IIB5p . 

First, we compare the most probable matching with the actual matching and analyze their dependence on the 
diffusivity n* . We define the Hamming distance between the actual and the most probable matchings as the fraction 
of edges which are present in the actual matching but are missed in the most probable matching. The comparison is 
shown in Fig. [51 where we plot the Hamming distance as a function of the diffusivity. The inset of the Figure shows a 
transition observed at k* < 0.174(4) from the lower diffusivity phase, where no difference between the most probable 
and actual matching were detected, to the high diffusivity phase where the most probable and actual matchings are 
distinctly different. The Hamming distance saturates to 1 as k* — > oo. Note that this phase transition is a property 
specific to the random distance model and that it was not observed in our particle simulations in 2d and 3d discussed 
in the main text. 

To estimate the most probable value of the diffusion constant k based on the observed data (set of mutual distances) 
we need to maximize the full likelihood Z{k), eq. (2), with respect to k, i.e. to minimize the free energy of the model. 
In the limit, N ^ oo, the minimum of the free energy within the random distance model is always achieved at k = k*. 
However, the minimum gets flatter as the the actual diffusivity k* grows. In a finite size system the statistical error 
in estimating k* is 1/ \J {NJ-"). Let us remind here that if the actual matching is known then the statistical error 
is At* (2/N). In Fig. [SI we plot (in log- log-scale) the second derivative of the free energy as a function of k* . We 
observe that for k* < 1 the second derivative is roughly T" = 1/[2(k*)^], hence in this regime the statistical error 
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FIG. 8: Left: Result of the zero temperature population dynamics with A'pop = 10000, c — 20, tmax ~ 10000. Distance between 
the actual and the most probable matchings shown as a function of diffusivity ft* . Inset: Zoom into the region of the phase 
transition at low diffusivity. At about Kc ~ 0.174(4) the most probable matching ceases to be equal (asymptotically) to the 
actual matching. Right: The second derivative of the free energy at k* as a function of k*. We observe the change of the slope 
at ft* « 1 from 1/[2(k*)^] to a larger value. 



is of a similar origin as if the actual matching would be known. For k* > I, the scaling of the second derivative 
changes to roughly proportional to 1/k*, hence in this regime estimation of the actual k* from the full likelihood 
is more accurate. It is noteworthy to emphasize that the crossover in the second derivative (and thus accuracy of 
the log- likelihood based prediction) is observed at k* w I, whereas the phase transition in the Hamming distance 
(between the actual and most probable matchings) takes place at the much lower values k* = Kc ~ 0.174(4). 



Appendix C: Comparative analysis of BP and MCMC algorithms 

In the body of the paper we stated the problem of learning the parameters of the flow/diffusion from particle 
tracking data in terms of the partition function of an associated matching problem. One important point articulated 
in the manuscript is that, even in the regime where the most probable matching is very different from the actual 
one, computing the partition function enables us to estimate the parameters accurately. Note, that to the best of 
our knowledge all the methods used so far in particle tracking reconstruction/learning, have been relying solely on a 
single matching. 

Our results show that BP is both efficient and accurate for approximating the maximum of the partition function. 
We have used an MCMC-learning algorithm, but our original focus was primarily on assessing the quality of the BP 
by comparing it to the MCMC near-to-cxact result. To guarantee near-to-exact result from the MCMC algorithm the 
running time is impractical. In this Section we give some additional details and comment also on comparative speed 
performance of the two learning algorithms. 

MCMC sampling is the most common method for estimating properties of the Boltzmann distribution. However, 
standard MCMC approaches sample from the given distribution aiming to evaluate certain expectation values and 
they are not directly suitable for accurate counting (evaluation of the partition function) required for our purposes. A 
considerable speed up in sampling of lower cost configurations is usually achieved via Metropolis-Hastings implemen- 
tation of the MCMC, however, in order to compute the partition function (or the free energy) one has to numerically 
integrate the energy function over a range of temperatures, and this integration is computationally costly. 

To resolve this problem we have used a special MCMC algorithm uniquely designed for this purpose - a somewhat 
simplified version of the Fully Polynomial Randomized Approximation Scheme (FPRAS), originally introduced for 
computing permanents (and this is what our partition function is) in [l3l . The unabridged version of this FPRAS 
algorithm is guaranteed to give a value that is not worse than (1 + e)L [L being the exact value) in 0{N^^) and 
poly(l/e) computational time [iTj . A description of our simplified algorithm follows. 

Let us denote the matrix whose permanent we want to compute by A. In our particular case, A contains the 
likelihoods of pairs of particles in different time steps corresponding to each other. The algorithm works in stages. It 
starts with a constant matrix A', whose elements are all equal to the largest value in A, amax, and whose per- 
manent is easily computed as {amax)^n\. In each stage i, elements of A' are reduced by a constant factor to 
= max{exp(— l/2)a^^, a„„}, thus bringing A' closer to A. Using S (a parameter) samples from all possible 
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perfect matchings (i.e., permutations), an approximate ratio rj is computed of the new value of the permanent of A' 
and its old value (before reducing elements of A' in step i). The samples are obtained using a Metropolis-Hastings 
MCMC, where neighbors of a state are all matchings with any two pairings swapped, and the probability of accepting 
a proposed step is the ratio between products of values of A' corresponding to new and old matchings. This MCMC 
is run for T (a second parameter) steps for each sample. Finally, the algorithm finishes when the matrix A' becomes 
sufficiently close to the original A. The approximate value of the permanent is then obtained as Yii^i ' io,max)"'n\. 
The two parameters, S and T, determine the time complexity of the algorithm, which is 0{ST). 
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FIG. 9: Logarithm of the partition function computed from the belief propagation algorithm (in red thick line) and from the 
Monte Carlo described in the text. The Monte Carlo is guaranteed to find an e approximation of the exact value in 0{N^^) 
time. We, however, run it for a much shorter periods and analyze how does its performance depend on the running time. We 
find that in order to compare with BP in the estimation of the maximum, MCMC needs simulations that are at least one order 
of magnitude longer. 



To facilitate comparison of the MCMC and BP efficiency (speed) , we downgraded strict quality requirement of the 
original FPRAS-MCMC and conducted simulations to access its deterioration in quality with the computational time 
(and number of samples) decrease. For this comparison we have used instances of the particle tracking described 
in Section [Xj The results, shown in Fig. |9l lead to conclusion that for comparable performance we need about an 
order of magnitude longer running time for the MCMC. Note, however, that we tested only specific implementation 
of the MCMC algorithm, and obviously further comparative analysis of BP-based, MCMC-based and possibly other 
learning schemes will be necessary. 

We find it useful to complete this Section by discussing another important point related to the comparison of the 
MCMC and BP efficiency for evaluating the maximum of the partition function. This is a crucial quantity for learning 
as it provides the most likely estimate of the flow/diffusion parameters. We have observed that using BP provides 
an explicit estimate not only for Z, but also for its first and second derivatives over the parameters. This observation 
has helped us to significantly accelerate the BP-based learning thus streamlining the search for the maximum of the 
partition function over the parameters. The acceleration is achieved via the use of Newton method in parallel with BP 
iterations. The trick itself gives an advantage to BP and it is currently not clear if similar accelerations arc possible 
to achieve in the MCMC-based learning. 



Appendix D: Possible extensions 



In this section we discuss possible extensions of our approach that will be interesting to pursue in future work. 
We nevertheless find it useful to briefly discuss them here, in order to give a better idea about possibilities (and 
limitations) of a systematic probabilistic approach to particle tracking. 

Incorporating uncertainty Our approach allows simple modiflcations to incorporate various types of uncer- 
tainty. For example, one can account for particles leaving and entering the box (focal volume in the experimental 
realization) via the following relaxation of the constraint C ({c}) in eq. (1), 



c(w)^n 




n 



(Dl) 
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where the parameter fj, is interpreted as a chemical potential that accounts for the fluctuations in the number of 
particles. The BP equations can be rewritten for this form of likelihood and the penalty parameter /i would have to 
be varied so as to optimize the likelihood and allow appropriate number of non-matching particles. 

Processing of multiple snapshots Our approach can also be adopted to inference and learning based on 
more than two subsequent snapshots. If the incremental displacement of particles between two subsequent snapshots 
depends only on particle positions at the early snapshot in the pair and does not depend on other details of the 
preceding evolution of the particles, then generalization of our approach consists in treating matchings between each 
pair of snapshot in the sequence independently, stating the total likelihood of the sequence as a product of the 
respective expressions eq. (1) for the pairs, and making global optimization over the governing parameters 6 for the 
entire product. 

If the motion of particles has memory spanning over m snapshots, then our approach can still be implemented, 
however this will require knowledge or estimation for the joint probability of observations P{xi, . . . , Xm\0). Note that 
in this correlated case already the subtask of inferring MPA assignment, equivalent to the so-called mult i- matching 
problem [l^, |3l|, becomes significantly more difficult than in the non-correlated case. 

Some sort of predictor-corrector scheme could also be investigated to solve the case with memory, or to achieve a 
speed-up of multiple snapshot processing. These promising possibilities are yet to be investigated and tested. 

Use of marginal probabilities We concentrated on inference of fiow parameters in regime where tracking of 
individual particles is no longer possible. Note, however, that once the correct parameters are estimated the belief 
propagation provides also estimates of probabilities with which a given particle from the first snapshot moved to a 
given position in the second snapshot. This information may be useful in cases where we are interested in information 
about the trajectories of individual particles. The values of these probabilities also provide concrete information about 
the level of uncertainty where a given particle moved. 

Note, however, that if a particle is matched to the position of its most probable image, the resulting configuration 
may even not be a one-to-one mapping. Indeed, a position in the second image might be the most probable one for 
multiple particles from the first image. Decimation-like techniques, of the type used to solve similar inconsistency 
problems in the constraint satisfaction problems [s^ . can then be used to obtain a consistent one-to-one mapping. In 
either case, even though the resulting configuration may not be equal to the actual displacement, it may still provide 
a useful visual information. 

Learning multi-scale flow In this paragraph we describe how the learning framework, illustrated in the body 
of the manuscript on the model case of the Batchelor (single-scale) velocity, can be extended to the more general 
case of a multi-scale flow. The Batchelor model of velocity considered in the main text assumed that s, the matrix 
of velocity gradients, does not depend on the particle index, i.e. in other words, all the particles in the cloud sense 
the same velocity gradient. In a more realistic turbulence setting the velocity gradient varies on the spatial scale 
correspondent to the viscous (Kolmogorov scale). Thus, when particles are seeded in a cloud exceeding in size the 
viscous scale, one needs to use a richer higher-parametric model, for example stated in terms of the following harmonic 
expansion of the instantaneous velocity field: u{r) = exp(iqr)wq, where the number of harmonics (number of 
terms in the sum) is Ng. Then, velocity advecting particle i, where labeling is according to the first of the two 
consecutive images, is modeled in terms of the following gradient fiow, Vi{ri\r*) = Si{xi){ri — r*) + vf ^r*), where 

vf\r*) = u{r*) = 6^P(*'3"'!*)"q' — V"u'^(ri) = X]q '^^P(*9^r)*'Z""q! ^^'^ ^j* i'' ^ pre-selected and 

time-independent (frozen) grid point closest to the particle's initial position (in the first image), Xi = ri{t = 0). (One 
possible choice for r* might be Xi itself.) The stochastic dynamic equation for trajectory of particle i becomes 

n = vf\rl) + sf^-'^^'ivDin - r*) + Ut). (D2) 

These equations can be integrated over time, thus resulting in the following generalization of Eqs. (7-9) from the main 
text 

P^ix,,y3) = (detAf,)-5exp(-if"(Af-i)fV) ; (D3) 
f = y^- - r* - W,(Xf dt'Wr\t)vl"^ + - r*); (D4) 
M, = nW,iA)[j^^W,'\t)Wr'-^it)dt] W^iA), (D5) 

where Wi{t) = exp(isi). The task of reconstruction becomes to infer Ng harmonics given positions of N particles 
in the two snapshots, for example assuming that the diffusion coefficient k is known. Naturally, a reliable multi- 
parametric reconstruction is feasible if Ng ^ A^. Extension of our BP-based scheme to this multi-parametric setting 
is straightforward, as the graphical model employed to solve the problem is identical to the one described in the main 
text for the Batchelor model case. The main technical difficulty in learning the multi-parametric velocity will in fact 
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be in maximizing the partition function of the model over a larger (than in the Batchelor case) number of the degrees 
of freedoms (harmonics). However, and as discussed above in Section [CJ this problem can be dealt with efficiently 
by alternating Newton (parameter adjustment) steps with BP iterative steps. The results and performance of the 
multi-scale inference are, however, yet to be tested. 

Interacting particles Let us note that extending our approach to systems of interacting particles may be possible, 
but it would constitute a more significant challenge, as in this case possible probabilistic models of particles' evolution 
arc much harder to evaluate. 
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